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ABSTRACT 

We consider the effects of eccentricity on the fragmentation of gravitationally unstable accretion disks, using 
numerical hydrodynamics. We find that eccentricity does not affect the overall stability of the disk against 
fragmentation, but significantly alters the manner in which such fragments accrete gas. Variable tidal forces 
around an eccentric orbit slow the accretion process, and suppress the formation of weakly-bound clumps. The 
"stellar" mass function resulting from the fragmentation of an eccentric disk is found to have a significantly 
higher characteristic mass than that from a corresponding circular disk. We discuss our results in terms of the 
disk(s) of massive stars at ~ O.lpc from the Galactic Center, and find that the fragmentation of an eccentric 
accretion disk, due to gravitational instability, is a viable mechanism for the formation of these systems. 
Subject headings: Galaxy: center - stars: luminosity function, mass function - accretion, accretion disks - 
hydrodynamics - methods: numerical 



1. INTRODUCTION 

The relative proximity of the Galactic Center (henceforth 
GC) provides a unique opportunity for "close-up" study of the 
processes that influence the formation of galaxies and black 
holes. Recent advances in telescope technology have enabled 
the resolution of individual stars in the crowded GC envi- 
ronme nt, and these new data have produced several puzzles 
(e.g., Ghez et al."1998t iGenzel et allllOOl iGhez et al.ll2005t 
[Pau mard et al. 2006). In particular, the detection of large 
numbers of young, massive stars is strongly suggestive of re- 
cent star formation activity at or close to the GC, despite the 
fact that the extreme physical conditions at the GC (temper- 
ature, tidal shear, etc.) prohibit the formation of stars by the 
same mechanisms that form stars in the solar neighborhood. 
The data show that a small number (few tens) of B-type stars 
exist very close to the GC (within ~ O.Olpc; these are the so- 
called "S-stars"), and a larger population (> 100) of massive 
O, B and Wolf-Rayet (WR) stars exist at slightly larger radii 
(~ 0. Ipc). Many of this second popu lation of stars are k nown 
to form a coherent ring or disk (e.g.. lGenzel et al.ll2003l) and, 
while the existenc e of a second such disk of stars rem ains 
controversial (e.g.. iPaumard et al.1 120061: iLu et al.ll2006l) . the 
existence of at least one stellar disk seems secure. 

A popular theory posits that the stellar disk(s) formed 
via fragmentation of an accretion disk due to grav- 
itational instability ( e.g.. | Polvachenko & ShukmanI Il977t 



Levin & Beloborodo 



)lJ2003f 

m 



Shlosman & Begelman 1987, 



1989t iGoodman & TanlEooS lNavakshinl i2006[ iLevinii2007i) . 
However, we have recently demonstrated that the ob- 
served eccentricities of the stellar orbits (typically > 0.3, 
iBeloborodov et al]l2006l) are inconsistent wi th dynamical re- 
laxati on from a circular initial configuration (lAlexander et al.l 
l2007h . unless the stellar mass function (MF) is extremely 
top-heavy, much more so than is inferred from o bservations 
(iNayakshin & Sunyaevll2005HPaumard et al.1 '2006). We sug- 
gested that fragmentation of an eccentric accretion disk could 
explain this discrepancy, and here seek to investigate the ef- 
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fects of eccentricity on disk fragmentation. 

In recent years gravitational instabilities in accretion disks 
have been the subject of much theoretical and numerical re- 
search, primarily in the context of prot oplanetary disks an d 
planet formation (see, e.g., the review by'Duri sen et al.l2007h . 
Some details remain contested, but the basic physical pro- 
cesses are now well understood. In order for a gravitationally 
unstable disk to fragment, two conditions must be satisfied. 
Firstl y, the disk m ust be sufficiently massive and/or cold that 
the iToomrd (1 1 964) criterion is satisfied: 



CsJ7 ^ 
^ ttGS ^ 



(1) 



Here c^ is the local sound speed, fl the local orbital frequency, 
and S the disk surface density. This criterion ensures that the 
disk is gravitationally unstable, but in order for the instabil- 
ity to lead to fragmentation the disk must also cool rapidly, 
or else compre ssional h eating will stabilize the disk against 
fragmentation. iGammid (|2001) used a local analysis to show 
that the cooling timescale of the disk must satisfy 



fcooi ^ 3ri 



(2) 



and subsequent work using global simulations has verified 
the validity of this criterion (Rice et al. 2003; Lodato & Rici 
|2004 l2005h . The criteri on also depends weakly on th e 
adopted equation of state (iGammid 120011: iRice et al.ll2005h . 
and studies of protoplanetary disks have shown that this de - 
pendence can be non-trivial in real disks (B olev et al]|2007h . 
In the case of an eccentric disk, the circularization energy pro- 
vides an additional, and potentially important, source of heat- 
ing, and this remains largely unstudied to date. 

Theoretical models of disk fragmentation around black 
holes typically consider circular disks with near-equilibrium 
initial conditions, implicitly assuming that such a config- 
uration results from some sort of "accre tion event" (such 
as the capture of a molecular clou d, e . g .. iNavakshinI 120061 : 
INayakshin et"ani2007t lLevinll2007h . In the case of the GC, 
the so-called circumnuclear ring at a few pc from Sgr A* is 
thought to be the likely gas reservoir for such accretion events 
(lMorrislll993t ISanderslll998h . but the dynamics of the infall 
process is not well understood. Moreover, some models of 
cloud infall predict the formation of coherent, eccentric disks 
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around the central black hole (e.g. JSandersll 19981) . We expect 
such disks to be gravitationally unstable, and star formation 
via fragmentation of an eccentric disk would provide an el- 
e gant solution to the pro blem of the stellar orbits discussed 
in lAlexander et al. ( 2007), although it may leave some issues, 
such as the z-distribution of the orbits, unresolved. In this 
paper, therefore, we seek to understand the dynamics of grav- 
itational instability and fragmentation in eccentric disks. In 
Section|2]we describe our numerical simulations, and present 
the results in Section [3] We discuss the implications of our 
results for models of star formation at the GC, as well as the 
limitations of our analysis, in Section HI and summarize our 
conclusions in Section|5] 

2. SIMULATIONS 

Our simulations make use of the publicly-available 
smoothed-parti cle hydrodynamics (SPH) code GADGET2 
(ISpringelll2005h . The code has been modified to allow the use 
of a simple cooling prescription (as outlined below), which 
is valid only in a disk geom etry. We adopt the standard 
iMonaghan & Gingoldl (Il983h prescription for the S PH vis- 
cosity (with aspH = 1 -0), using the "Balsara-switch" (iBalsaral 
Il995h to limit the artificial shear viscosity (as specified in 
Equations 11-12 of Springel 2005), and adopt sufficiently 
high numerical resolution that the transport of angular mo- 
mentum due to numerical dissipation is much less than that 
expected from self-gravitating angular m omentum transport 
jLodato & Rice' l2004l2005tlNelsonll200 6). As demanded by 
iNelson (2006), we allow for a variable gravitational soften- 
ing length, and fix the SPH smoothing and gravitational soft- 
ening lengths to be equal throughout. We use the standard 
Barnes-Hut for malism to compute the gravitational force tree 
(as described in ISpringell2005 ). and use A^neb = 64 ± 2 as the 
number of SPH neighbours. Lastly, in order to avoid close 
orbits around the "black hole" limiting the time-step to an 
unreasonably small value, we use a single sink particle as 
the central gravitating mass. This particle simply accretes all 
gas particles that p ass within its sink radius (as described in 
ICuadra et al.ll2006h . and is merely a numerical convenience 
that has no physical effect on the simulation. We set the sink 
radius of the "black hole" to be 1/4 of the inner disk radius (or 
semi-major axis, if e ^Q), which has a value of 1 in our scale- 
free simulations (see Section l23T l. In practice fewer than 0.1% 
of the particles are swallowed by the sink particle in any of 
our simulations. The simulations were run on the TUNGSTEN 
Xeon Linux cluster at NCSA^, using 64 parallel CPUs for the 
highest-resolution runs. 

2.1. Initial conditions 

We set up our initial disk as follows, u sing an approach that 
closely mirrors that of iRice et al.l (l2005h . We adopt a surface 
density profile 



E(a) oc a 
and a disk temperature that scales as 
T{a) oc fl~ ' , 



(3) 



(4) 



where a is the orbital semi-major axis. The disk sound speed 
C.5 oc r'/^ is normalized so that the Toomre parameter 



Q- 



c,n 



ttGS 
' See |http: //www.ncsa.uiuc. edu/| 



(5) 



is equal to 2.0 at the outer boundary. (Here il is the orbital fre- 
quency.) Using this prescription, the Toomre parameter scales 
approximately as 

Q oc a-^l^ (6) 

in the initial disk. With this set-up the initial disk is marginally 
gravitationally stable, and cools into instability with time. 

In order to set up such a disk in SPH, we first define a 
disk mass Mj, and set the mass of each gas particle to be 
m = Md/Ns,pH, where A^sph is the number of gas particles. We 
then divide the disk into concentric annuli, and distribute the 
particles in the annuli according to the surface density pro- 
file. Within each annulus, particles are distributed such that 
the mass flux around the annulus is constant. For an orbit of 
given semi-major axis a, eccentricity e, and azimuthal angle 
(p, the radius is given by 



K0) = 



1- 



l+ecos( 



(7) 



and the radial and azimuthal components of the velocity are 
given by 



GM,, 



Vr = e sm < 



V0: 



r(l+ecos(j)) 

GMenc(l+eCOS(/)) 



(8) 



(9) 



respectively, where Menc is the mass enclosed by the orbit 
(the mass of the central object plus the mass of the disk at 
semi-major axis < a). Note that, as our disks are sufficiently 
massive to be self-gravitating, the disk mass can cause a sig- 
nificant perturbation to the Keplerian potential of the central 
mass. (Strictly these expressions only apply to a spherically 
symmetric mass distribution, but they are sufficiently accurate 
for our purposes.) The particles are then distributed vertically 
by randomly sampling a hydrostatic equilibrium (Gaussian) 
density distribution with scale-height H = c^/Uir) (although 
we note that this is not strictly an equilibrium configuration 
whene ^0). 

2.2. Cooling 

In order to draw comparisons with previous studies, we 
adopt a simple cooling prescription of the form 

(10) 



cool./ 



Ui 
^cool 



where m, is the internal energy of particle /, and the cooling 
timescale fcooi is given by 



t -^ 



(11) 



Here the constant /3 is an input parameter: previous studies 
have found that the fragmentation boundary in circular disks 
lies at (3 c± 3, with a weak dependerice on the eq uation of state 
adopted ( lGammiell200Tt iRice et aljlloolllool . 

In the case of an eccentric disk, however, there is an am- 
biguity in this definition of the cooling timescale. One can 
choose either to use fl(r) = y/GMuJr^, or n,(a) = y/GM\,h/a^. 
The first case results in a cooling timescale which varies 
around the orbit, while in the second case fcooi is constant for 
individual orbital streamlines. It is not entirely clear which of 
these is most appropriate when compared to a real disk, as this 
depends on the coolants involved. However, it seems unlikely 
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Simulation 


A'SPH 


e 


/3 


^max 


^^dZ-'Wbh 


Mill, smoothing 
length? 


Wruns 


CIRC 


500,000 





3 


25.0 


0.2 


No 


3 


CIRCB5 


500,000 





5 


25.0 


0.2 


No 




ECC.25 


500,000 


0.25 


3 


25.0 


0.2 


No 




ECC.50 


500,000 


0.5 


3 


25.0 


0.2 


No 




ECC.50B5 


500,000 


0.5 


5 


25.0 


0.2 


No 




ECC.50VAR 500,000 


0.5 


3 


25.0 


0.2 


No 




ECC.75 


500,000 


0.75 


3 


25.0 


0.2 


No 




CIRCHR 


5x 10« 





3 


5.0 


0.05 


Yes 




ECCHR 


5x 10* 


0.5 


3 


5.0 


0.05 


Yes 





TABLE 1 
List of simulations run, showing resolution and physical pai'ameters for each 
simulation. A'sph is the number of SPH particles used to model the gas disk, 

and Nium is the number of different random realizations of the initial 
conditions that were used in each case. The simulation ECC.50VAR used the 

cooling prescription that results in a cooling time which varies around the 

eccentric orbit; all the others used a cooling time that was constant along the 

orbital streamlines. 

that the rate at which an individual fluid element cools will 
vary on a timescale shorter than the dynamical timescale, so 
we adopt the second prescription for most of our modelling. 
We do, however, run a model with varying cooling time as a 
test case. 
Consequently, we define the cooling timescale as 



^cool ~ 



n(a) 



= (3\ 



GMkh' 



(12) 



and compute the semi-major axis a of a given particle directly 
from the instantaneous orbital elements (based on the assump- 
tion of Keplerian orbits). We adopt an adiabatic equation of 
state throughout, with adiabatic index 7 = 5/3. 

2.3. Simulation details 

In order to study the effects of eccentricity on disk fragmen- 
tation, we have conducted a number of simulations. The pa- 
rameters of these simulations are specified in Table [1] and are 
summarized below. The simulations are scale-free: we adopt 
a system of units where the length unit is the semi-major axis 
of the inner disk edge, the mass unit is that of the central black 
hole, and the time unit is the orbital period at a = 1 . 

Firstly, we conducted a large suite of simulations at moder- 
ate resolution, using 500,000 gas particles. In order to resolve 
the fragmentation process properly we adopted a rather mas- 
sive (and therefore thick) disk, with Md/Mbh = 0.2, and simply 
stopped the calculations when the density of collapsed regions 
became high enough to restrict the timestep to an unfeasibly 
small value. We adopted a dynamic range of 25 in semi-major 
axis for these simulations. For a circular disk (e = 0) we con- 
ducted three identical simulations with (3 = 3 (i.e., that are 
marginally unstable to fragmentation), each using a different 
random realization of the initial distribution of gas particles 
(i.e.,, a different noise field), and similarly performed three 
such simulations with e = 0.5. In addition, we performed a 
single simulation with e = 0.25, another with e = 0.75, and a 
further simulation with e = 0.5 using the cooling prescription 
that varies around the disk orbit (see Section l2!2] i. Lastly, we 
performed two simulations with a longer cooling time of (3 = 5 
(whic h should be stab le against fragmentation in the circular 
case; Rice et al.ll2003 [): one with e = 0, and one with e = 0.5. 

As seen in Section |372l below, these low -resolution simula- 
tions allow us to infer a great deal about the effects of eccen- 
tricity on disk fragmentation, but do not run for long enough, 
or produce enough fragments, to allow measurement of a MF 



Consequently, we also ran two further simulations at much 
higher resolution, using 5 x 10^ gas particles and a somewhat 
smaller radial range (with the outer edge of the disk at 5 times 
the semi-major axis of the inner edge). We adopted the light- 
est disk that could be resolved safely throughout (at a factor 
of ~ 2 better than the Nelson 2006 criterion for spatial resolu- 
tion in the vertical direction), and so adopted Mj/Mbh = 0.05. 
If we scale this to the GC system, this gives an SPH particle 
mass of O.O3M0 (assuming a black hole mass of 3 x lO^'M©). 
We note that this is slightly more massiv e than the disks 
adopt ed by previous such simulations (e.g., [Nayakshin et ajj 
120071 who adopt Md/Mbh = 0.01), and also slightly more mas- 
sive than current observational estimates of the total stellar 
mas s in the GC disk(s) (wh i ch range from c ± 0. 003-0. 03Mbh, 
e.g.. iNavakshin et al.ll2006t iPaumard et al.l l2006). However, 
such observations may under-estimate the mass of the ini- 
tial gas disk as, give n the age of the stellar disk(s) (~ 6Myr, 
iPaumard et al.ll2006l) . it is reasonable to assume that the most 
massive stars that formed are no longer present. Moreover, 
given the scale-free nature of our simulations, it is unlikely 
that this small over-estimate of the disk mass will have a sig- 
nificant effect on our results. 

Additionally, in order to follow the fragmentation process 
beyond the initial collapse phase in these high-resolution sim- 
ulations, we imposed a minimum SPH smoothing length of 
0.001 . This is well below the Hill radius expected for even the 
least massive clumps, and merely acts to limit collapse to den- 
sities that would stop the calculation; the clumps essentially 
have density profiles which are "flat-topped" within this ra- 
dius. We prefer this to a sink particle approach because bound 
material can become unbound from fragments as they move 
around an eccentric orbit (due to the variable shear), and sink 
particl es of the type used by, for example, INavakshin et al.l 
(120071) cannot lose mass in this manner (by construction). We 
conducted two of these high-resolution simulations: one of a 
circular disk, and one with e = 0.5. 

3. RESULTS 

3.1. Code Tests 

As a test of our numerical method, we first set out to re- 
produce previously obtained results. In the circular case, 
we find that the simulation with /3 = 5 was indeed stable 
against fragmentation, and instead produced quasi-steady an- 
gular momentum transport via spiral density waves. This sim- 
ulation was allowed to run for 5 outer orbital periods (i.e., 625 
inner orbital periods) and showed no evidence of fragmenta- 
tion, suggesting that this configuration is indeed stable against 
fragmentation over long timescales. By contrast, all of the 
simulations with (3 = 3 were unstable to fragmentation, typ- 
ically breaking into fragments after around 90 inner orbital 
periods'*. Reflections of density waves from the inner edge 
of the disk act to stabilize the region within a few H of the 
boundary, so that the first fragments form at radii of ~ 3-5 
(where the orbital timescale is approximately 10 times that at 
the inner boundary). Thus we find that our numerical proce- 
dure was successful in reproducing previous results, namely 
that the fragmentation boundary for a gravit ationally unsta- 
ble disk with 7 = 5/3 lies at fcooi^ — 3 (e.g., lGammiell2001l; 

■* The exact time at which fragmentation occurs is somewhat arbitrary, as 
it merely reflects the time required for the initial configuration to cool into 
instability (the disk initially has Q ~ 20 in the inner regions). However, the 
fact that different random realizations of the initial conditions result in frag- 
mentation at the same time supports our assertion that the simulations are 
numerically converged. 
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iRice et aljllool . 

3.2. Low-resolution runs 

The low resolution runs make use of a rather massive (and 
therefore thick) disk and, if they fragment at all, are not ex- 
pected to produce very many fragments (as the most unstable 
length scale is always of order H). Consequently, the anal- 
ysis of these low resolution runs is limited to if, when and 
where they fragment, and the number of fragments produced: 
there are too few fragments to allow detailed analysis of their 
properties. In these runs the density is allowed to increase 
without limit, so the simulations simply stop when the den- 
sity becomes sufficiently high to force a very small timestep. 
Typically this occurs 1-2 (local) orbital periods after the first 
fragments form, and only a small fraction of the gas is ac- 
creted into clumps before the simulations stop. 

In order to compare the fragmentation properties of runs 
with different eccentricities, we require a simple method to 
quantify the behavior (masses, velocities, positions, etc.) of 
the fragments that form. We define fragments (or "clumps") 
to be bound objects containing at least 128 SPH particles 
(i.e., with at least double the number of nearest neighbors). 
The clumps are first identified by locating peaks in the den- 
sity distribution, and are then tested for bound-ness. For sim- 
plicity we calculate the potential energy of pairs of particles 
using a slightly simplified gravitational potential of the form 
Gm/{r+h), where h is the local SPH smoothing length: this is 
sufficiently accurate for our purposes, and is faster to compute 
than convolution with the rather complex smoothing function 
used by GADGET2. The masses of individual clumps are de- 
termined by ranking particles by total energy (potential plus 
kinetic plus thermal) and iterating outwards until the first un- 
bound particle is reached. We define the position and velocity 
of each clump as the mean values of all of the SPH parti- 
cles bound to that clump. While somewhat cumbersome, we 
prefer this method of identifying clumps to a geometric one 
(using, for example, concentric shells) as it allows accurate 
treatment of non-spherical clumps, and we apply this method 
identically to all of our simulations (both here and in the high- 
resolution runs discussed below). 

3.2.1. Fragmentation boundary for e =/0 

Our first aim was to study the effect of eccentricity on the 
overall stability of the disk against fragmentation. The en- 
ergy liberated by circularization of an eccentric orbit (with a 
fixed angular momentum) is given simply by the product of 
e^ and the total energy of the corresponding circular orbit, so 
in principle this represents a large energy reservoir that could 
stabilize the disk against fragmentation. However, in our sim- 
ulations this energy is not liberated sufficiently quickly to 
prevent fragmentation: all of the eccentric disks with /? = 3 
fragmented (apart from that with e = 0.75, see below), while 
none of the models with (3 = 5 showed any evidence for frag- 
mentation despite running for several outer orbital times (see 
FiglH'^- Moreover, the fact that the disk with e = 0.75 and 
P = 3 did not fragment is likely an artefact of our numeri- 
cal set-up. In this simulation the major- to-minor axis ratio 
of the disk is so large that when the inner disk precesses it 
approaches the outer edge of the disk. At this point density 
waves are reflected back off the outer edge of the disk, and 

^ Visualiz ations of our SPH simulations were created using SPLASH: see 
iPncg J2007h for details. 
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Fig. 1 . — Snapshots of disk midplane density in the central region of the 
low-resolution runs with e = 0.5: the upper panel shows a model with /3 = 
3 (simulation ECC.50), the lower /3 = 5 (simulation ECC.50B5). Both are 
plotted on the same (logarithmic) density scale, and both are shown for t = 
112.5. The /3 = 3 run has fragmented, while the /3 = 5 run remains stable. The 
major axis of the disk was initially aligned with the ,t-axis: the precession of 
the inner disk is clearly seen in both snapshots. 

appear to stabilize the disk against fragmentation. In prin- 
ciple this could be a physical effect, but it is not clear that 
a real disk would have such a sharp outer edge with such a 
large major-to-minor axis ratio. Additional simulations with 
a much larger dynamical range in semi-major axis (and there- 
fore a much larger particle number) are required to investigate 
this further. We therefore conclude that the disk eccentric- 
ity has little or no effect on the location of the fragmentation 
boundary, except possibly at very high eccentricity (e > 0.75). 
In a sense, this can be attributed to our choice of initial con- 
ditions. Our disks were set up so that the angle of pericenter 
was constant with radius, and consequently none of the ini- 
tial streamlines intersect. We expect the circularization en- 
ergy to be liberated by shocks, which arise when differential 
precession of the disk causes different orbital streamlines to 
intersect. In our models the central body provides a point- 
mass potential, and the only non-point-mass contribution to 
the potential (i.e., the part that causes precession) comes from 
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'Low resolution' simulations: N = 250,000 , Mj/M^j,=0.2 



Evolution of clump mass 




t (after first fragment) / inner orbital period 

Fig. 2. — Evolution of the number of bound clumps in simulations with 
different eccentricity: the time axis is normalised to the point where the first 
bound clump forms in each simulation. The solid black (e = 0) and red (e = 
0.5) fines sfiow tfie mean of thi'ee different simulations; the (noisier) dashed 
green line (e = 0.25) is for a single simulation only. The eccentric disks 
clearly form fewer bound clumps than are seen in the circular case. The 
decfine in clump number at f ~ 15 in the circular disks is due to mergers. 

the disk itself. The precession timescale is approximately 
(■^bh/-^d)forb^ but the disk becomes unstable to fragmentation 
after approximately 1-2 local cooling times. Consequently 
our disks become unstable to fragmentation before they are 
able to liberate a significant fraction of their circularization 
energy (as seen in FiglT), and the eccentricity has no signifi- 
cant effect on the overall stability criteria. We find that a disk 
with uniform eccentricity is just as likely to be unstable to 
fragmentation as a circular disk with the same thermodynamic 
properties, but note that our disks are artificially constructed 
to have such a uniform eccentricity. The question of whether 
or not these initial conditions are realistic is discussed in Sec- 
tion |4lT] 



3.2.2. The effects of disk eccentricity on the fragmentation 

process 

While it seems that eccentricity has little or no effect on the 
overall stability of a self-gravitating disk, our simulations do 
show that the details of the fragmentation process are strongly 
affected by eccentricity. This is clearly seen in Fig|2] which 
shows the number of bound clumps as a function of time for 
simulations with different eccentricity. Although all of the 
disks with /3 = 3 fragment, clump formation is clearly sup- 
pressed in the cases where e ^0, suggesting that disk eccen- 
tricity has a significant influence on the accretion of gas on 
to individual clumps. This is presumably due to the variable 
tidal forces around an eccentric orbit: gas which is bound to a 
particular clump at apocenter is not necessarily bound at peri- 
center In the simulations with e y^Q, the most weakly-bound 
clumps do not survive their first pericenter passage, and in- 
stead are sheared apart by tidal forces. The simulation with a 
variable cooling time (simulation ECC.50VAR, not shown in 
Fig|2|i shows less pronounced tidal stripping than those with 
constant cooling time (as the cooling is faster near pericen- 
ter), but still results in many fewer fragments than seen in the 
circular disks. 

The effect of tidal forces is highlighted further in Fig|3] 
which shows the mass of a particular clump in one of the 
e = 0.5 runs during the first orbit following its formation. The 




Orbital phase / 2n 

Fig. 3. — Evolution of the mass of a typical clump, formed in a disk with 
e = 0.5 (simulation ECC.50). The horizontal axis shows the orbital phase, 
with the zero-point taken at pericenter; the vertical axis shows the mass bound 
to the clump, in units of the mass of the central black hole. The plot begins at 
the point where the clump first becomes bound, and ends where the simula- 
tion stops (giving a total period of approximately one orbit for this particular 
clump). The decrease in mass caused by tidal stripping during pericenter pas- 
sage is clearly visible. The small variations in mass close to apocenter arise 
when the clump passes through local enhancements in the gas density, caused 
by global spiral density waves. 

mass of the clump is approximately constant around much of 
the orbit, but decreases by around 70% as the clump passes 
pericenter This clearly demonstrates that the variable tidal 
forces in an eccentric disk have a strong effect on the accre- 
tion process, and that tidal stripping of this kind must be taken 
into account when modelling this process. 

3.3. High-resolution runs 

In order to study the effect of eccentricity on the frag- 
mentation process in more detail, we conducted two addi- 
tional simulations wit h mu ch higher numerical resolution. As 
described in Section 12.31 these simulations used 5 million 
SPH particles and Md/M^h = 0.05, resulting in an SPH par- 
ticle mass of 10~^Mbh. In addition, we im posed a minimurn 
SPH smoothing length on these simulations dMonaghanl 1 992t 
iNavakshin et al.ll2007 ), limiting collapse at the highest densi- 
ties and allowing us to follow the simulations for much longer 
after the formation of the first fragments. The thinner disk 
used here results in many more fragments being formed, and 
following the simulations for longer allows us to construct 
MFs of the formed clumps with good number statisti cs. We 
note that the recent work of iPickett & DurisenI (l2007h found, 
using an isothermal equation of state, that the artificial vis- 
cosity used in SPH can increase the survival times of weakly- 
bound clumps in simulations such as these. In principle this 
could influence our results, so we adopt a resolution limit 
which ensures that all our clumps are well-resolved, and this 
should minimize the effects of artificial dissipation on our de- 
rived mass functions. We use the algorithm described in Sec- 
tion l3.2l above to identify clumps, but reject clumps with fewer 
than 192 bound SPH particles (i.e., 3A^ngb) from our subse- 
quent analysis as such clumps are only marginally resolved. 
Scaled to the GC system, this corresponds to a resolution limit 
of ~ 6M0 in mass and ~ 20AU in length (assuming a black 
hole mass of 3 x 1O*M0 and a characteristic disk radius of 
O.lpc). We compute the number of bound clumps for every 
output snapshot (every 0.05 time units), but for computational 
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Fig. 4. — Final MF for the circular disk, measured at the point where the 
simulation stopped (t = 12.7). The solid histogram shows the measured MF; 
the dashed lines show the Salpeter slope and best-fitting log-normal function 
(with a width of 0.79dex). The vertical dotted fine denotes the resolution 
limit. 

reasons we only evaluate the MF for a handful of these snap- 
shots (every 1 time unit, starting from the point where the first 
bound clump forms). 

3.3.1. Circular disk 

The reference simulation of a circular disk (CIRCHR), ran 
for 12.7 inner orbital periods. The first bound fragment was 
formed at f = 6.2, and most of the fragmentation occurred at 
radii < 2 (simply because the dynamical time in the outer disk 
is too long for fragmentation to occur). At the point where 
the calculation stopped it had formed 313 bound clumps, of 
which 263 were above the (mass) resolution limit. Within 
r = 2 the bound clumps account for 40.5% of the total disk 
mass, suggesting that much of the gas has been accreted into 
bound objects and the resulting MF is indeed representative. 
The final MF is shown in FiglD The slope is close to Salpeter 
at high masses, but the MF shows a w ell-resolved turnover a t 
~ 10"^Mbh and is in fact well-fit by a lMiller & ScaJol (Il979h 
log-normal distribution function, with a width of ~ O.Sdex 
and a characteristic mass of 1.1 x 10~^Mbh- Note that this 
is essentially the Toomre mass, which for this simulation is 
ttYjH^ ci 10~^Mbh: this corresponds to a characteristic mass 
of ~ 3OM0 when scaled to the GC. The MF is clearly not 
consistent with the Salpeter slope at even moderate stellar 
masses, and the low-m ass turnover is consisten t with previous 
numerical simulations jNavakshin et al.ll2007h. and also with 
observations of the GC ste llar disks (iNayakshin & SunyaevI 
l2005t iPaumard et alJl2006l) . Thus we are confident that our 
reference model is both numerically accurate and physically 
plausible, and now look to the differences that result when the 
disk is eccentric. 

3.3.2. Eccentric disk 

The simulation ECCHR ran for 13.5 inner orbital periods, 
and the first fragment formed at f = 6.75. As in the circular 
disk, the fragmentation was confined to a < 2, and both sim- 
ulations ran for ~ 7 inner orbital periods after the formation 
of the first fragment. As we saw in the low resolution runs, 
the eccentric disk formed fewer clumps: at the point where 
this calculation stopped it had formed 189 clumps, of which 
178 had masses above the resolution limit. Within a = 2 the 




Fig. 5. — Snapshots of disk midplane density in the central region of the 
high-resolution run with e = 0.5 (simulation ECCHR). The snapshot is taken 
at the point where the simulation ended, t = 13.5, and many bound clumps 
are clearly visible. The major axis of the disk was again initially aligned with 
the .v-axis: the less massive disk results in considerably less precession than 
was seen in the low-resolution runs. 

bound clumps had accreted 26.5% of the available gas. The 
clumps follow orbits very similar to the disk: while there is 
some scatter in the instantaneous values of the clump eccen- 
tricities, the rms eccentricity of all of the bound clumps at 
the end of the simulation was ~ 0.49 (essentially equal to the 
disk eccentricity of 0.5). The final midplane density profile is 
shown in Fig|5] 

3.3.3. Mass function evolution 

In order to look at the effects of eccentricity on the frag- 
mentation and accretion processes, it is instructive to look 
at the evolution of the MF in both the circular and eccentric 
disks. Fig|6] shows the evolution of both clump number and 
accreted mass in the inner part (a < 2) of both the circular 
and eccentric disks. Both disks show a "two-phase" forma- 
tion process: a fragmentation phase, where the clump num- 
ber increases rapidly and most of the bound mass is in small 
clumps; followed by an accretion phase, where the clump 
number remains approximately constant but the clumps con- 
tinue to accrete gas from their surroundings. Once again we 
see that the eccentric disk forms fewer clumps than the cir- 
cular disk, and we again attribute this to the tidal destruction 
of weakly-bound clumps in the eccentric disk. Moreover, the 
accretion of mass is much slower in the eccentric disk, sug- 
gesting that tidal forces affect the accretion process across the 
full range of clump masses. This is reinforced by the mea- 
sured MFs: the eccentric disk MF consistently shows a larger 
characteristic clump mass than seen the circular disk. 

The MFs are compared directly in Fig|7] measured at a 
point where both simulations have accreted the same fraction 
of gas into bound clumps. The eccentric disk MF is taken 
at the end of the simulation (when 26.5% of the disk mass 
has been accreted into clumps), and the "reference" circular 
disk MF is taken 4.4 time units after the formation of the first 
clump (26.8% accreted). Both are reasonably well-fit by log- 
normal distribution functions: the best-fitting log-normal in 
the circular case has a characteristic mass of 6.6 x 10~* and a 
width of 0.6 dex; in the eccentric case the best-fitting charac- 
teristic mass is 1 .2 x 10"^, and the width is 0.5 dex. A KS test 
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Fig. 6. — Evolution of bound clumps in the high-resolution simulations. 
The solid lines show the evolution of the number of bound clumps (left-hand 
scale), while the points show the fraction of disk mass within a < 2 accreted 
into bound clumps (right-hand scale). The dotted lines between the points are 
added as a guide to the eye. As in Fig|2] the time axis is normalized to the 
point where the first bound clump forms in each simulation. 

gives a probability of 7 x 10~* that the two MFs are drawn 
from the same underlying distribution, and similar compar- 
isons at earlier times in the simulations (with the same ac- 
creted mass fractions) all result in probabilities of < 1 % that 
the two MFs are the same. Thus we conclude that the disk 
eccentricity alters the accretion process significantly, primar- 
ily because of tidal effects around the eccentric orbit. Mass is 
tidally stripped from bound clumps during pericenter passage, 
with two notable effects. Firstly, the most weakly-bound, 
lowest-mass clumps fail to survive pericenter passage, result- 
ing in a relative dearth of low-mass objects in the clump MF 
Secondly, the overall accretion rate is slowed, suggesting that 
an eccentric disk will take significantly longer than a circular 
one to accrete all of its mass into bound clumps. 

4. DISCUSSION 

While our simulations are scale-free, they have obvi- 
ous applications to the formation of the stellar disk(s) ob- 
served at the GC. Observations of the GC stellar disks find 
that the stellar MF must be much more "top-heavy" than 
the standard Salpeter MF. The relatively low X-ray flux 
from the GC region suggests a significant deficit of low- 
mass (< 5M(7)) stars compared to local star-forming regions 
(JNavakshin & Sunvaev 2005), and estimates of the stellar MF 
from the /T-band luminosity function suggests a MF slope that 
is 1.0 -1.5 dex flatter than the Salpeter slope (iPaumard et alJ 
|2006|) . Also, given th at the total stellar mass in the disk(s) 
is_~ IO'^-IO^Mq fe.g..'Navak shin et al.ll2006l: iPaumard et al.' 
l2006l) . the large number (> 100) of observed massive O- 
& WR-stars suggests a characteristic stellar mass of ~ 10- 
50Mq, much larger than that observed in the solar neigh- 
bourhood. This is consi stent with our results ( and those of 
previous studies such as Nay akshin et al.ll2007h . which pre- 
dict a characteristic stellar mass of ~ SOM© (when scaled 
to the GC), and a MF that is significantly deficient in low- 
mass objects. Previous studies have not been able to repro- 
duce all of the observed properties of the disks (notably the 
eccentric stellar orbits), but our results suggest that fragmen- 
tation of moderately eccentric accretion disk can form disks 
of stars with masses and orbits consistent with those observed. 
Our simulations suggest that eccentric disks fragment to pro- 
duce MFs that are slightly more top-heavy than expected in 
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Fig. 7. — Comparison of clump MFs in the circular (black) and eccentric 
(red) disks. As in Fig|4] the dotted line shows the resolution limit, and the 
dashed line the Salpeter slope. The MF for the eccentric disk is measured at 
the end of the simulation; in the circular case it is measured where the same 
fraction of the disk mass has been accreted into bound clumps. The eccentric 
disk clearly forms fewer low-mass objects than the circular disk, and a KS test 
rejects the possibility that the two MFs are drawn from the same underlying 
distribution. 



the circular case (due to the tidal destruction of the lowest- 
mass clumps), but not so top-heavy as to be inconsistent with 
the observed stellar populations. Moreover, previous simula- 
tions have shown that changes in the thermodynamic proper- 
ties of the disk, such as the cooling rate, can alter the MF at 
a level at least as signi ficant as that seen in our simulations 
dNavakshin et al.l 120071) . Consequently, we are satisfied that 
our results are robust, within the range of parameter space ex- 
plored by our simulations. However, there are three major 
areas in which our simulations differ from real systems, and 
we discuss each of these in turn below. 

In addition, we note that angular momentum transport (by 
spiral waves) is negligible in the simulations that exhibit frag- 
mentation, as the fragmentation and subsequent accretion of 
gas occurs on the dynamical timescale (which is much shorter 
than the timescale for angular momentum transport). It has 
long been thought that angular momentum transport due to 
gravitational instabilities can be responsible for significant 
accretion in disks around black holes (e.g., Shlosma n et alJ 
19 89). This is indeed li kely, but our results confirm those 
of iNavakshin et al] (l2007i) . that star formation due to gravita- 
tional instability in accretion disks occurs on a timescale that 
is significantly shorter than the timescale for angular momen- 
tum transport. Consequently, if star formation occurs at all 
radii in the disk, it seems likely that little or no gas will be ac- 
creted on to the central black hole. We note, however, that real 
disks around black holes are like ly only unstable to fragmen- 
tation in their outer regions (e.g.. lLevinll2007l;lKing & Pringld 
l2007h . and that whether or not gravitational instability in such 
disks results in fragmentation or transport depends critically 
on the issues of cooling discussed in Section l4~2l 

Our simulations are most applicable to the GC stellar 
disk(s), but their scale-free nature allows us to consider the re- 
sults in terms of other systems also. It has been suggested that 
the young stellar population at the GC may be characteristic 
of more generi c episod es of black hole (AGN) activity i n other 
galax ies (e.g., iLevin & Beloborodov..2003j: Thompson et al.l 
l2005h . and this has impUcations for models of AGN accre- 
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tion and the growth of supermassive black holes. The conse- 
quences of our models for AG N fe eding depend on the disk 
thermodynamics (see Section 14.2b . but we note in passing 
that our models fit naturally into, for example, the black hole 
growth scenario presented by King & Pringle (2006, 2007]). 
The applicability of our results to protoplanetary discs, how- 
ever, is less clear It seems likely that protoplanetary disks 
are g ravitationally unstable at early times in their evolution 
(e.g.. iDurisen et al] |2007). but it is less clear that such disks 
possess significant eccentricities. If eccentric, self-gravitating 
disks do exist around young stars then our results could easily 
be applied to such systems, but at present there is no strong 
evidence for the existence of such systems. 

4.1. Initial configuration 

The first area of concern is in our initial conditions. We 
choose an initial configuration that is marginally stable, and 
allow the disks to cool into instabilit y. This ap proach has been 
followed by previous studies (e.g., iRice et al. 2003, 2005), 
and recent work suggests that use of these artificial initial 
conditions does not affect the f ragmentation boundary in cir- 
cular disks (IClarke et al.ll2007l) . What is less clear, however, 
is whether or not our choice of a uniform angles of pericen- 
ter is similarly valid. As discussed in Section [T] in systems 
like the GC we expect the disk to form via some ty pe of ac- 
cretio n event, and while previous simulations (e.g., ISandersI 
Il998l) have found that an eccentric configuration can result, 
there is no particular reason to expect a uniform eccentricity. 
If the disk does not have uniform eccentricity then neighbor- 
ing orbital streamlines will intersect, and the resulting shocks 
have the potential to provide significant additional heating. 
Such heating may act to stabilize the disk against fragmen- 
tation, but a realistic treatment of this problem requires full 
(magneto-)hydrodynamic simulations of the infall and cap- 
ture of a gas cloud by the black hole, and such modeling is 
clearly beyond the scope of this work. What we have demon- 
strated, however, is that if an eccentric disk is able to evolve 
towards a configuration where it satisfies the (circular) condi- 
tions for fragmentation, then it will fragment before any such 
shock heating is able to stabiUze the disk. 

4.2. Cooling and thermodynamics 

The second major area of uncertainty in our models is 
our simplified coo ling prescription. Our cooling prescrip- 
tion (Section I2.2l l parametrizes the cooling rate by semi- 
major axis only, and while this is ideal for the sort of nu- 
merical experiments conducted here, it is far from realis- 
tic. In the case of protoplanetary disks, the subject of heat- 
ing and cooling rates in gravitationally unstab le disks is an 
area of on-going debate (e.g., the review by iDurisen et alj 
l2007h . but in the case of a disk around a black hole 
we can justify the approximate form of our cooling pre- 
scription by a few simple scaling arguments. Similar ar- 
gume nts have been discussed previously by other authors 
(e.g.. Levin & Beloborodov"2003^, 'Johnson & Gammie||2003t 
iGoodman & Tan 2004; Levin 2007; Rafikov 2007), but we re- 
state them here for completeness. 

A self-gravitating disk has Q^l, and therefore if we know 
the disk surface density we can estimate the sound speed, and 
therefore the temperature, in the disk. If we scale our model 
to the GC system (by adopting a characteristic radius of 0. Ipc 
and Mbh = 3 x lO^'Mo), we find a typical disk surface den- 
sity of ~ lOOg cm"^ and a corresponding disk temperature 
of ^ lOOK (estimated at a = 2, but these quantities do not 



vary significantly over the limited radial range of our simu- 
lations). At such temperatures the dominant source of opac- 
ity is dust grains, and typical val ues for the Ross eland mean 
opacity are k ~ 1-lOcm^g"' (e.g.. lSemenov et al] r2003). The 
opacity does not change dramatically at temperatures below 
that at which the dust grains sublime (typically 2± 1500K). 
As seen above, scaling our models to the GC system (assum- 
ing that the initial gas mass in the disk is comparable to the 
observed stellar mass) results in a much lower temperature 
than this, so significant variations in the disk opacity (such as 
the opacity gap discus sed by, e.g..l Johnson & Gammiell2003l; 
[Thompson et al.ll2005h seem unlikely. Only if star formation 
in such disks is very inefficient, contrary to what is seen in our 
(and other) simulations, would we expect the disk tempera- 
tures to be high enough for dust sublimation to be significant. 
Therefore, we expect such a disk to be marginally optically 
thick^, and the resulting cooling timescale is ty pically compa- 
rable to the orbital period (e.g.. lRafikovll200'fl) . However, we 
also note that our cooling model takes no account of changes 
in the cooling function as the local density increases (due to 
gravitational collapse). A more realistic model would make 
use of an opacity-based cooling model that accounts for the 
density structure of the disc, and would also adopt a somewhat 
stiffer equation of state in the the collapsing clumps. The ap- 
plication of such models to protoplanetary disks is still mired 
in controversy (e.g., Boley et al. 2006; Mayer et al. 2007), but 
the physical conditions in a GC disk are known to be much 
more conducive to fragmentati on than those in protoplanetary 
disks (see e.g., the discussion in lRafikovl2007l) . Consequently 
we expect that our disks will cool rapidly enough to be unsta- 
ble to fragmentation. Moreover, the uncertainties discussed 
here are not specific to our models of eccentric disks, so we 
do not expect our results regarding the differences between 
eccentric and circular disks to be affected by these simplifica- 
tions. 

A further simplification is that our cooling prescription de- 
pends only on the local conditions in the disk, and is not in- 
fluenced at all by the global disk structure. Previous work 
using this scale-free, local cooling prescription has found that 
it gives rise to instabilities that are well-approximated by a 
local description (Lodato & Rice 2004, 2005), and therefore 
the global disk structure is unlikely to have a significant ef- 
fect on the spatial distribution of clumps in our simulations. 
However, in a real disk the cooling time is not scale-free and, 
in general, gravitational instabilities cannot be well-described 
locally (Balbus & Papaloizou 1999). Numerical simulations 
of protoplanetary disks have found that transport can be dom- 
inated by low-order global rnodes i n the non-scale-free case 
(iBolev et al.1 120061: ICai et all IIoOTL IPurisen et al.1 lIOOTh . It 
seems likely, however, that the local approximation will be 
more accurate for the very thin disks considered here than in 
the case of protoplanetary disks (where M^/M^. '^ 0.1). 

What is less certain is whether or not a realistic system is 
able to evolve into the configuration discussed here. The virial 
temperature of a gas cloud that falls towards the GC from, 
for example, the circumnuclear ring is orders of mag nitude 
larger than the disk temperatures considered here (e.g Morris] 
1993; Sanders 1998) . An infalling cloud may not virialize 
before it forms a coherent disk structure, but there is clearly 
sufficient energy available to heat the gas significantly. At 
temperatures from ~ 1000-10"^K the typical Rosseland opac- 

* Note that this also verifies our earlier assumption that the cooling rate is 
unlikely to change significantly around the orbit. 
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ities are much smaller than those at ~ lOOK, as the primary 
coolants (dust, and also molecular speci es) are destroyed at 
such temperatures (i Semenov et al.ll2()03h . Thus, we expect 
the cooling rates at higher temperatures to be much lower than 
those considered here, and it may be difficult for such a disk 
to cool into i nstability in a more realisti c model (see also the 
discussion in lJohnson & Gammiell2003h . Our simple thermo- 
dynamic model is viable once the disk reaches the state con- 
sidered here but, as discussed above, we do not address the 
issue of whether or not a real system will be able to reach 
such a configuration. 

4.3. Accretion physics at stellar scales 

The last major area neglected by our simulations is the 
small-scale physics of star-formation. Our simulations have 
sufficient numerical resolution to treat the fragmentation pro- 
cess correctly, but when scaled to the GC our minimum length 
resolution is a few tens of AU. This is 2-3 orders or magnitude 
larger than the typical stellar radius, and we make no attempt 
to model the physics of star formation on small scales. Es- 
sentially our MFs are MPs of pre-stellar cores, and as such 
they may not be representative of the stellar MFs that result 
from disk fragmentation. However, we note that the fragmen- 
tation into bound clumps occurs on orbital timescales, and the 
results presented in Fig|6] suggest that all of the disk gas will 
be bound to clumps within ~ 10 orbital periods. The orbital 
period of the GC stellar disks is ^ lOOOyr, but the typical 
timescales for such cores to collapse to form stars are ^ 10"*- 
lO^yr (e.g., McKee & Tan 2003). Consequently we do not ex- 
pect stellar feedback processes to have a significant effect on 
the process of disk fragmentation. However, individual bound 
clumps may well form more than one star eac h, and we als o 
expect to see some mergers between clumps (lLevinll2007h . 
More detailed simulations, using higher numerical resolution 



and more realistic thermal physics, wiU be required to address 
these issues further. 



5. SUMMARY 

We have constructed hydrodynamic models of gravitation- 
ally unstable disks, and investigated the effects of eccentricity 
on disk fragmentation. We find that the fragmentation bound- 
ary is unaffected by eccentricity, and that conditions which 
lead to fragmentation in circular disks do likewise in eccen- 
tric disks. However, we find that the fragmentation process is 
altered by eccentricity, as the variable tidal forces around an 
eccentric orbit have a strong effect on the accretion of gas on 
to bound fragments. We find that the formation of low-mass 
fragments is suppressed, and that the growth (via accretion) 
of more massive fragments proceeds much more slowly in an 
eccentric disk than in a circular one. We consider our results 
in the context of the GC stellar disks, and find good agreement 
between the MFs in our simulations and those observed. We 
find that fragmentation of eccentric accretion disks is a viable 
mechanism for the formation of the GC stellar disks, with the 
resulting stellar MF and orbits consistent with observations. 
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